//============================================================================
// Name        : FourierLaw.cpp
// Author      : Carlos Olivares
// Version     :
// Copyright   : My copyright is Love
// Description : Hello World in C, Ansi-style
//============================================================================

//Its simple just need a chain, with potentials and noise all with the same mass 
// First the deterministic version just a chain to see how they move
// Even more basic just a harmonic body
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <time.h>

double fAnharmonic(double k1, double k2, double x){
	return -k1*x - k2*x*x*x;
}

int main(int argc, char* argv[]) {
   // FIXME: Input arguments by command line.
   if(argc != 8){
		printf("Usage: FourierLaw [prefix_file] [N_chain] [k1] [k2] [total_time] [dt] [T1] [T2] [lambda]\n");
		return EXIT_FAILURE;
	}
   
   //----------- Input Parameters.
   char *prefix_file;
   prefix_file	=	argv[1];
	const int N	=	atoi(argv[2]);
   double k1 = atof(argv[3]);
   double k2 = atof(argv[4]);
   int ttime_total = atoi(argv[5]);
   double dt =	atof(argv[6]); //0.01;
   double lambda = atof(argv[7]); 
   if (lambda <= 0) {
		printf("lambda must be positive\n");
		return EXIT_FAILURE;
	}
   double mu=dt/(lambda); // mu = 1/lambda
   // TODO: Add temperature effects -- dahh, of course!!!
//   double T1 = atof(argv[5]);
//   double T2 = atof(argv[6]);
   
   
   //----------- Output data file.
   char filename_pattern[90];
   char filename_dat[100];
   char ext_dat[]=".dat";
   
   strcpy(filename_pattern,prefix_file);
	for (int i=2;i<argc;i++){
		strcat(filename_pattern,"_");
		strcat(filename_pattern,argv[i]);
	}
   strcpy(filename_dat,filename_pattern);
   strcat(filename_dat,ext_dat);
   //char filename_dat[]="WithoutShotNoise.dat";
   
   printf("%d\t%f\t%f\n",ttime_total,dt,k1);
   FILE *fp_dat;
   fp_dat = fopen(filename_dat,"w");
   fprintf(fp_dat,"#$t$\t$x$\t$v$\trandom\n");
   
   
   //----------- Initial conditions.
	double x_old=1, v_old=0;
	double x_new=0, v_new=0;
	double ttime=0;
   
   //----------- Random number generation.
   const gsl_rng_type *T; T = gsl_rng_default;
	gsl_rng *r; r = gsl_rng_alloc(T);
	long unsigned int seed;
	seed = time(NULL);
	gsl_rng_set(r,seed);
   
   
   int l=0;
   int delta_l=0;
   int sign_il =1;
   
   printf("Molecular dynamics begin\n");

   //----------- Molecular dynamics.
	for (int i = 0; i < ttime_total; i++){
      
      
      // Poisson random change of velocity direction
      if (l==i) {
         
         sign_il = -1;
         delta_l=gsl_ran_poisson(r,mu);
         l = l+delta_l;
      }
      else{ 
         sign_il = 1;
      }
      
      // Motions Equations
		v_new = sign_il*v_old  + fAnharmonic(k1, k2, x_old)*dt;
		x_new = x_old + v_new*dt;

		x_old = x_new;
		v_old = v_new;
		ttime = i*dt;
		fprintf(fp_dat,"%f\t%f\t%f\t%d\n",ttime,x_new,v_new,sign_il);
	}
   
   fclose(fp_dat);
	return EXIT_SUCCESS;
}
